# Reproduction files - Study 2 descriptives, opinion formation

# Deepening Bridging and Moving Minds in Stressful Times







# Session Information -----------------------------------------------------




#R version 4.2.2 (2022-10-31 ucrt)
#Platform: x86_64-w64-mingw32/x64 (64-bit)
#Running under: Windows 10 x64 (build 19045)
#
#Matrix products: default
#
#locale:
#  [1] LC_COLLATE=German_Germany.utf8  LC_CTYPE=German_Germany.utf8    LC_MONETARY=German_Germany.utf8 LC_NUMERIC=C                   
#[5] LC_TIME=German_Germany.utf8    
#
#attached base packages:
#  [1] stats     graphics  grDevices utils     datasets  methods   base     
#
#other attached packages:
#  [1] car_3.1-1       carData_3.0-5   patchwork_1.3.0 scales_1.3.0    ggeffects_1.1.5 stargazer_5.2.3 nnet_7.3-18     sjPlot_2.8.12  
#[9] psych_2.2.9     forcats_1.0.0   stringr_1.5.0   dplyr_1.1.4     purrr_1.0.1     readr_2.1.4     tidyr_1.3.0     tibble_3.2.1   
#[17] ggplot2_3.5.1   tidyverse_1.3.2 haven_2.5.3     sjmisc_2.8.9   
#
#loaded via a namespace (and not attached):
#  [1] nlme_3.1-160        fs_1.6.1            lubridate_1.9.2     insight_0.19.0      httr_1.4.4          tools_4.2.2         backports_1.4.1    
#[8] R6_2.5.1            sjlabelled_1.2.0    DBI_1.1.3           colorspace_2.1-0    withr_2.5.0         tidyselect_1.2.0    Exact_3.2          
#[15] mnormt_2.1.1        emmeans_1.8.4-1     compiler_4.2.2      performance_0.10.2  cli_3.6.0           rvest_1.0.3         expm_0.999-7       
#[22] xml2_1.3.3          sandwich_3.0-2      bayestestR_0.13.0   mvtnorm_1.1-3       proxy_0.4-27        minqa_1.2.5         pkgconfig_2.0.3    
#[29] lme4_1.1-31         dbplyr_2.3.0        rlang_1.1.4         readxl_1.4.2        rstudioapi_0.14     farver_2.1.1        generics_0.1.3     
#[36] zoo_1.8-11          jsonlite_1.8.4      googlesheets4_1.0.1 magrittr_2.0.3      Matrix_1.5-1        Rcpp_1.0.10         DescTools_0.99.50  
#[43] munsell_0.5.0       abind_1.4-5         lifecycle_1.0.3     stringi_1.7.12      multcomp_1.4-22     MASS_7.3-58.1       rootSolve_1.8.2.3  
#[50] grid_4.2.2          parallel_4.2.2      crayon_1.5.2        lmom_2.9            lattice_0.20-45     splines_4.2.2       sjstats_0.18.2     
#[57] hms_1.1.2           knitr_1.42          pillar_1.10.1       ggpubr_0.6.0        boot_1.3-28         gld_2.6.6           estimability_1.4.1 
#[64] ggsignif_0.6.4      codetools_0.2-18    reprex_2.0.2        glue_1.6.2          data.table_1.14.6   modelr_0.1.10       nloptr_2.0.3       
#[71] vctrs_0.6.5         tzdb_0.3.0          cellranger_1.1.0    gtable_0.3.1        assertthat_0.2.1    datawizard_0.6.5    xfun_0.37          
#[78] xtable_1.8-4        broom_1.0.3         e1071_1.7-13        coda_0.19-4         rstatix_0.7.2       survival_3.4-0      class_7.3-20       
#[85] googledrive_2.0.0   gargle_1.3.0        timechange_0.2.0    TH.data_1.1-1       ellipsis_0.3.2 




# Calculations ------------------------------------------------------------







# load necessary packages
library(sjmisc)
library(haven)
library(tidyverse)
library(psych)
library(sjPlot)
library(nnet)
library(stargazer)
library(ggeffects)
library(scales)
library(patchwork)
library(car)



# 1. Data preperation, descriptives ---------------------------------------

# read dataset
data_AT <- readRDS("data_at_wght.RDS")


# number of cases
nrow(data_AT) 
frq(data_AT$finished_post) # 2140 completes in austria

data_AT %>% 
  filter(!is.na(position_post_cont & !is.na(position_pre_cont))) %>% nrow()

# descriptive information pre and post survey position
describe(data_AT$position_pre_cont) # 6.43
describe(data_AT$position_post_cont) # 6.54

# create new variable: pro health vs. pro liberties
data_AT$pos_pre_group <- car::recode(data_AT$position_pre_cont,
                                     recodes = "1:5 = 0; 6:10 = 1", as.factor = T)
frq(data_AT$pos_pre_group)
# 37.46 pro freedom; 62.54 pro health

# new variable strong prior opinion
data_AT$strong_prior <- car::recode(data_AT$position_pre_cont,
                                    recodes = "c(1,2,9,10) = 1; 3:8 = 0", as.factor = T)
frq(data_AT$strong_prior)
# 69.46 % strong prior; 30.54 % no strong prior


# diversity of the sample

# age
describe(data_AT$age) # mean = 44.5

# gender
frq(data_AT$gender) # 58.8 % female, 40.9 % male



# 2. Opinion change -------------------------------------------------------

# descriptive information pre and post survey opinion
describe(data_AT$position_pre_cont) # mean pre: 6.43
describe(data_AT$position_post_cont) # mean post: 6.54

# aggregate opinion change
data_AT <- data_AT %>%
  mutate(abs_position_change = abs(data_AT$position_post_cont - data_AT$position_pre_cont),
         position_change1 = ifelse(abs_position_change <= 0, "no change", "change"),
         position_change2 = ifelse(abs_position_change <=1, "no change", "change"))

# average individual opinion shift
mean(data_AT$abs_position_change, na.rm = T) # avg. change 0.29 scale points


# calculate absolute diff
data_AT$abs_change = abs(data_AT$position_post_cont - data_AT$position_pre_cont)

# mean
mean_abs_change <- mean(data_AT$abs_change, na.rm = TRUE)

# sd
sd_abs_change <- sd(data_AT$abs_change, na.rm = TRUE)

# Cohen's d 
cohen_d_abs <- mean_abs_change / sd_abs_change
cohen_d_abs







# individual changes
frq(data_AT$position_change1) # 16.18 % changed by at least 1 SP
frq(data_AT$position_change2) # 5.44  % changed by at least 2 SP


# direction of opinion change
data_AT <- data_AT %>%
  mutate(diff_pre_post = position_post_cont - position_pre_cont,
         opinion_formation = ifelse(
           abs_position_change == 0, "Stability", ifelse(
             diff_pre_post > 0 & position_pre_cont>=6, "Polarization", ifelse(
               diff_pre_post < 0 & position_pre_cont>=6, "Depolarization", ifelse(
                 diff_pre_post > 0 & position_pre_cont <= 5, "Depolarization", ifelse(
                   diff_pre_post < 0 & position_pre_cont <= 5, "Polarization", NA))))),
         opinion_formation = as.factor(opinion_formation))

frq(data_AT$opinion_formation) # 8.68 % depolarized; 7.5 polarized; 83.82 stability



# 3. Randomization check --------------------------------------------------

# check wether random split into 5 treatment groups worked


frq(data_AT$group_treatment) # cases are equally distributed among 5 groups

# age
ggplot(data_AT, aes(age, group_treatment, group = group_treatment))+
  geom_violin(draw_quantiles = 0.5)+
  geom_jitter(alpha = 0.1)
describeBy(data_AT$age, group = data_AT$group_treatment)
# visual: no differences

# paired t-test
TukeyHSD(aov(age ~ as.factor(group_treatment), data = data_AT))
# no differences between 5 groups regarding age


# education
sjt.xtab(data_AT$education, data_AT$group_treatment, show.col.prc = T)
# no differences between 5 groups regarding education


# pre survey position
ggplot(data_AT, aes(position_pre_cont, group_treatment, group = group_treatment))+
  geom_violin(draw_quantiles = 0.5)+
  geom_jitter(alpha = 0.1)
describeBy(data_AT$position_pre_cont, group = data_AT$group_treatment)
# no differences 

# paired t-test
TukeyHSD(aov(position_pre_cont ~ as.factor(group_treatment), data = data_AT))
# no differences between 5 groups regarding position_pre_cont


# knowledge
ggplot(data_AT, aes(index_corona_knowledge, group_treatment, group = group_treatment))+
  geom_violin(draw_quantiles = 0.5)+
  geom_jitter(alpha = 0.1)
describeBy(data_AT$index_corona_knowledge, group = data_AT$group_treatment)
# visual: no differences

# paired t-test
TukeyHSD(aov(index_corona_knowledge ~ as.factor(group_treatment), data = data_AT))
# no differences between 5 groups regarding knowledge


# leftright positioning
ggplot(data_AT, aes(leftright, group_treatment, group = group_treatment))+
  geom_violin(draw_quantiles = 0.5)+
  geom_jitter(alpha = 0.1)
describeBy(data_AT$leftright, group = data_AT$group_treatment)
# visual: no differences

# paired t-test
TukeyHSD(aov(leftright ~ as.factor(group_treatment), data = data_AT))
# no differences between 5 groups regarding leftright positioning


# 4. Regression models ----------------------------------------------------


# Study 2 Model A ---------------------------------------------------------


# set opinion stability as reference category
levels(data_AT$opinion_formation)
data_AT$opinion_formation <- relevel(data_AT$opinion_formation, ref = "Stability")

# levels of treatment group
data_AT$group_treatment <- as.factor(data_AT$group_treatment)
levels(data_AT$group_treatment) <- c("Contestatory", "Collaborative", "Open communication", "Information only", "No information (placebo)")
frq(data_AT$group_treatment)
# reorder factor
data_AT$group_treatment <- factor(data_AT$group_treatment, levels = c("No information (placebo)", "Information only", "Contestatory", "Collaborative", "Open communication"))


# missing values
frq(data_AT$opinion_formation) # 2127 valid opinion change; 8 NA
frq(data_AT$group_treatment) # 2135 valid group_treatment; 0 NA
frq(data_AT$strong_prior) # 2133  valid strong_prior; 2 NA
frq(data_AT$index_corona_knowledge) # 2135 valid opinion change; 0 NA
frq(data_AT$leftright) # 2133  valid leftright; 2 NA
frq(data_AT$trust_gov) # 2118  valid trust_gov; 17 NA
frq(data_AT$trust_exp) # 2126  valid trust_exp; 9 NA
frq(data_AT$index_cognition) # 2134  valid index_cognition; 1 NA
frq(data_AT$index_evaluation) # 2133   valid index_evaluation; 2 NA
frq(data_AT$index_accuracy) # 2135   valid index_accuracy; 0 NA


# descriptives continuous independent variables
describe(data_AT$index_corona_knowledge)
describe(data_AT$leftright)
describe(data_AT$trust_gov)
describe(data_AT$trust_exp)
describe(data_AT$index_cognition)
describe(data_AT$index_evaluation)
describe(data_AT$index_accuracy)

# listwise deletion
data_AT <- data_AT %>% 
  filter(!is.na(opinion_formation) & !is.na(group_treatment) & !is.na(strong_prior) &
           !is.na(index_corona_knowledge) & !is.na(leftright) & !is.na(trust_gov) &
           !is.na(trust_exp) & !is.na(index_cognition) & !is.na(index_evaluation) &
           !is.na(index_accuracy))

nrow(data_AT) # 2110 complete cases


# estimate model
at_pol <- multinom(opinion_formation ~ group_treatment + strong_prior + index_corona_knowledge +
                     leftright + trust_gov + trust_exp + index_cognition +
                     index_evaluation + index_accuracy, data = data_AT) 

at_pol_N <- nrow(at_pol$residuals) 
at_pol_N # 2110 complete cases

# coef plot study 2 model A
coef_plot_study2_mA <- plot_model(at_pol, type = "est", show.values = T, transform = NULL,
                                  auto.label = F, vline.color = "black", 
                                  value.size = 2.5, dot.size = 1.5, line.size = .5, 
                                  value.offset =.35, order.terms = c(1:24),
                                  axis.labels= c("Accuracy motivation", "Neeed for Evaluation",
                                                 "Need for Cognition", "Trust experts",
                                                 "Trust government", "Leftright",
                                                 "Corona knowledge", "Strong prior",
                                                 "Open Communication", "Collaborative",
                                                 "Contestatory","Information Only"))+
  labs(title = "(A) All Participants", caption = paste0("N = ", at_pol_N))+
  theme_sjplot()+
  scale_y_continuous(limits = c(-3,3))
coef_plot_study2_mA[["layers"]][[1]][["aes_params"]]$size <- 0.3
coef_plot_study2_mA[["layers"]][[1]][["aes_params"]]$alpha <- 0.5

# predicted probabilities
pred_at_pol <- ggemmeans(at_pol, terms = "group_treatment", typical = "mean", condition = c(strong_prior =0, gender_male = 0))
pred_at_pol$x <- factor(pred_at_pol$x, levels = c("Open communication", "Collaborative", "Contestatory", "Information only", "No information (placebo)"))

prob_plot_study2_mA <- pred_at_pol %>%
  ggplot(aes(predicted, x))+
  geom_point()+
  geom_segment(aes(x=conf.low, y =x, xend=conf.high, yend=x))+
  facet_grid(~response.level, scales = "free")+
  xlab("")+
  ylab("")+
  theme_bw()+
  scale_x_continuous(n.breaks = 4, labels = label_percent(accuracy = 1))+
  theme(axis.text = element_text(size =7))


study2_plotA <- coef_plot_study2_mA / prob_plot_study2_mA





# Study 2 Model B ---------------------------------------------------------


# filter dataset, keep only pre survey pro health respondents
frq(data_AT$pos_pre_group) # N = 1319

data_AT_health <- data_AT %>% 
  filter(pos_pre_group == 1)

nrow(data_AT_health) # 1319 pro health complete cases


# estimate model
at_pol_health <- multinom(opinion_formation ~ group_treatment + strong_prior + index_corona_knowledge +
                            leftright + trust_gov + trust_exp + index_cognition +
                            index_evaluation + index_accuracy, data = data_AT_health) 

at_pol_health_N <- nrow(at_pol_health$residuals) 
at_pol_health_N # 1319 complete cases

# coef plot study 1 model B
coef_plot_study2_mB <- plot_model(at_pol_health, type = "est", show.values = T, transform = NULL,
                                  auto.label = F, vline.color = "black", 
                                  value.size = 2.5, dot.size = 1.5, line.size = .5, 
                                  value.offset =.35, order.terms = c(1:24),
                                  axis.labels= c("Accuracy motivation", "Neeed for Evaluation",
                                                 "Need for Cognition", "Trust experts",
                                                 "Trust government", "Leftright",
                                                 "Corona knowledge", "Strong prior",
                                                 "Open Communication", "Collaborative",
                                                 "Contestatory","Information Only"))+
  labs(title = "(B) Pro Health", caption = paste0("N = ", at_pol_health_N))+
  theme_sjplot()+
  scale_y_continuous(limits = c(-3,3))+
  theme(axis.text.y = element_blank())
coef_plot_study2_mB[["layers"]][[1]][["aes_params"]]$size <- 0.3
coef_plot_study2_mB[["layers"]][[1]][["aes_params"]]$alpha <- 0.5

# predicted probabilities
pred_at_pol_health <- ggemmeans(at_pol_health, terms = "group_treatment", typical = "mean", condition = c(strong_prior =0, gender_male = 0))
pred_at_pol_health$x <- factor(pred_at_pol_health$x, levels = c("Open communication", "Collaborative", "Contestatory", "Information only", "No information (placebo)"))

prob_plot_study2_mB <- pred_at_pol_health %>%
  ggplot(aes(predicted, x))+
  geom_point()+
  geom_segment(aes(x=conf.low, y =x, xend=conf.high, yend=x))+
  facet_grid(~response.level, scales = "free")+
  xlab("")+
  ylab("")+
  theme_bw()+
  scale_x_continuous(n.breaks = 4, labels = label_percent(accuracy = 1))+
  theme(axis.text.y = element_blank())


study2_plotB <- coef_plot_study2_mB / prob_plot_study2_mB


# Study 2 Model C ---------------------------------------------------------


# filter dataset, keep only pre survey pro freedom respondents
frq(data_AT$pos_pre_group) # N = 791

data_AT_free <- data_AT %>% 
  filter(pos_pre_group == 0)

nrow(data_AT_free) # 791 pro freedom complete cases


# estimate model
at_pol_free <- multinom(opinion_formation ~ group_treatment + strong_prior + index_corona_knowledge +
                            leftright + trust_gov + trust_exp + index_cognition +
                            index_evaluation + index_accuracy, data = data_AT_free) 

at_pol_free_N <- nrow(at_pol_free$residuals) 
at_pol_free_N # 791 complete cases

# coef plot study 1 model C
coef_plot_study2_mC <- plot_model(at_pol_free, type = "est", show.values = T, transform = NULL,
                                  auto.label = F, vline.color = "black", 
                                  value.size = 2.5, dot.size = 1.5, line.size = .5, 
                                  value.offset =.35, order.terms = c(1:24),
                                  axis.labels= c("Accuracy motivation", "Neeed for Evaluation",
                                                 "Need for Cognition", "Trust experts",
                                                 "Trust government", "Leftright",
                                                 "Corona knowledge", "Strong prior",
                                                 "Open Communication", "Collaborative",
                                                 "Contestatory","Information Only"))+
  labs(title = "(C) Pro Civil-Liberties", caption = paste0("N = ", at_pol_free_N))+
  theme_sjplot()+
  scale_y_continuous(limits = c(-3,3))+
  theme(axis.text.y = element_blank())
coef_plot_study2_mC[["layers"]][[1]][["aes_params"]]$size <- 0.3
coef_plot_study2_mC[["layers"]][[1]][["aes_params"]]$alpha <- 0.5

# predicted probabilities
pred_at_pol_free <- ggemmeans(at_pol_free, terms = "group_treatment", typical = "mean", condition = c(strong_prior =0, gender_male = 0))
pred_at_pol_free$x <- factor(pred_at_pol_free$x, levels = c("Open communication", "Collaborative", "Contestatory", "Information only", "No information (placebo)"))

prob_plot_study2_mC <- pred_at_pol_free %>%
  ggplot(aes(predicted, x))+
  geom_point()+
  geom_segment(aes(x=conf.low, y =x, xend=conf.high, yend=x))+
  facet_grid(~response.level, scales = "free")+
  xlab("")+
  ylab("")+
  theme_bw()+
  scale_x_continuous(n.breaks = 4, labels = label_percent(accuracy = 1))+
  theme(axis.text.y = element_blank())


study2_plotC <- coef_plot_study2_mC / prob_plot_study2_mC


# Fig. complete
(coef_plot_study2_mA + coef_plot_study2_mB + coef_plot_study2_mC) /
  (prob_plot_study2_mA + prob_plot_study2_mB + prob_plot_study2_mC)
ggsave("study2_regression_opinion change_logit.png", height = 7, width = 12, dpi = 300)




# 5. Model Appendix -------------------------------------------------------


# estimate model A
at_pol_append <- multinom(opinion_formation ~ group_treatment, data = data_AT) 

at_pol_append_N <- nrow(at_pol_append$residuals) 
at_pol_append_N # 2110 complete cases

# coef plot study 2 model A
coef_plot_study2_mA_a <- plot_model(at_pol_append, type = "est", show.values = T, transform = NULL,
                                    auto.label = F, vline.color = "black", 
                                    value.size = 2.5, dot.size = 1.5, line.size = .5, 
                                    value.offset =.35,
                                    axis.labels= c("Open Communication", "Collaborative",
                                                   "Contestatory","Information Only"))+
  labs(title = "(A) All Participants", caption = paste0("N = ", at_pol_append_N))+
  theme_sjplot()+
  scale_y_continuous(limits = c(-3,3))
coef_plot_study2_mA_a[["layers"]][[1]][["aes_params"]]$size <- 0.3
coef_plot_study2_mA_a[["layers"]][[1]][["aes_params"]]$alpha <- 0.5

# predicted probabilities
pred_at_pol_append <- ggemmeans(at_pol_append, terms = "group_treatment", typical = "mean", condition = c(strong_prior =0, gender_male = 0))
pred_at_pol_append$x <- factor(pred_at_pol_append$x, levels = c("Open communication", "Collaborative", "Contestatory", "Information only", "No information (placebo)"))

prob_plot_study2_mAa <- pred_at_pol_append %>%
  ggplot(aes(predicted, x))+
  geom_point()+
  geom_segment(aes(x=conf.low, y =x, xend=conf.high, yend=x))+
  facet_grid(~response.level, scales = "free")+
  xlab("")+
  ylab("")+
  theme_bw()+
  scale_x_continuous(n.breaks = 4, labels = label_percent(accuracy = 1))+
  theme(axis.text = element_text(size =7))


study2_plotA_append <- coef_plot_study2_mA_a / prob_plot_study2_mAa


# estimate model B
at_pol_health_append <- multinom(opinion_formation ~ group_treatment, data = data_AT_health) 

at_pol_health_N_append <- nrow(at_pol_health_append$residuals) 
at_pol_health_N_append # 1340 complete cases

# coef plot study 2 model B
coef_plot_study2_mB_a <- plot_model(at_pol_health_append, type = "est", show.values = T, transform = NULL,
                                    auto.label = F, vline.color = "black", 
                                    value.size = 2.5, dot.size = 1.5, line.size = .5, 
                                    value.offset =.35, col = "#377EB8",
                                    axis.labels= c("Open Communication", "Collaborative",
                                                   "Contestatory","Information Only"))+
  labs(title = "(B) Pro Health", caption = paste0("N = ", at_pol_health_N_append))+
  theme_sjplot()+
  scale_y_continuous(limits = c(-3,3))+
  theme(axis.text.y = element_blank())
coef_plot_study2_mB_a[["layers"]][[1]][["aes_params"]]$size <- 0.3
coef_plot_study2_mB_a[["layers"]][[1]][["aes_params"]]$alpha <- 0.5

# predicted probabilities
pred_at_pol_health_append <- ggemmeans(at_pol_health_append, terms = "group_treatment", typical = "mean", condition = c(strong_prior =0, gender_male = 0))
pred_at_pol_health_append$x <- factor(pred_at_pol_health_append$x, levels = c("Open communication", "Collaborative", "Contestatory", "Information only", "No information (placebo)"))

prob_plot_study2_mBa <- pred_at_pol_health_append %>%
  ggplot(aes(predicted, x))+
  geom_point()+
  geom_segment(aes(x=conf.low, y =x, xend=conf.high, yend=x))+
  facet_grid(~response.level, scales = "free")+
  xlab("")+
  ylab("")+
  theme_bw()+
  scale_x_continuous(n.breaks = 4, labels = label_percent(accuracy = 1))+
  theme(axis.text.y = element_blank())


study2_plotB_append <- coef_plot_study2_mB_a / prob_plot_study2_mBa



# estimate model C
at_pol_free_append <- multinom(opinion_formation ~ group_treatment, data = data_AT_free) 

at_pol_free_N_append <- nrow(at_pol_free_append$residuals) 
at_pol_free_N_append # 791 complete cases

# coef plot study 2 model C
coef_plot_study2_mC_a <- plot_model(at_pol_free_append, type = "est", show.values = T, transform = NULL,
                                    auto.label = F, vline.color = "black", 
                                    value.size = 2.5, dot.size = 1.5, line.size = .5, 
                                    value.offset =.35,
                                    axis.labels= c("Open Communication", "Collaborative",
                                                   "Contestatory","Information Only"))+
  labs(title = "(C) Pro Civil-Liberties", caption = paste0("N = ", at_pol_free_N_append))+
  theme_sjplot()+
  scale_y_continuous(limits = c(-3,3))+
  theme(axis.text.y = element_blank())
coef_plot_study2_mC_a[["layers"]][[1]][["aes_params"]]$size <- 0.3
coef_plot_study2_mC_a[["layers"]][[1]][["aes_params"]]$alpha <- 0.5

# predicted probabilities
pred_at_pol_free_append <- ggemmeans(at_pol_free_append, terms = "group_treatment", typical = "mean", condition = c(strong_prior =0, gender_male = 0))
pred_at_pol_free_append$x <- factor(pred_at_pol_free_append$x, levels = c("Open communication", "Collaborative", "Contestatory", "Information only", "No information (placebo)"))

prob_plot_study2_mCa <- pred_at_pol_free_append %>%
  ggplot(aes(predicted, x))+
  geom_point()+
  geom_segment(aes(x=conf.low, y =x, xend=conf.high, yend=x))+
  facet_grid(~response.level, scales = "free")+
  xlab("")+
  ylab("")+
  theme_bw()+
  scale_x_continuous(n.breaks = 4, labels = label_percent(accuracy = 1))+
  theme(axis.text.y = element_blank())


study2_plotC_append <- coef_plot_study2_mC_a / prob_plot_study2_mCa



# Fig. appendix complete
(coef_plot_study2_mA_a + coef_plot_study2_mB_a + coef_plot_study2_mC_a) /
  (prob_plot_study2_mAa + prob_plot_study2_mBa + prob_plot_study2_mCa)
ggsave("study2_regression_opinion change_appendix_logit.png", height = 7, width = 12, dpi = 300)



# descriptives study 2
frq(data_AT$strong_prior)
describe(data_AT$index_corona_knowledge)
describe(data_AT$leftright)
describe(data_AT$trust_gov)
describe(data_AT$trust_exp)
describe(data_AT$index_cognition)
describe(data_AT$index_evaluation)
describe(data_AT$index_accuracy)







# 1) Bivariate relationship between communication formats and justification rationality
data_AT_treat <- data_AT %>% filter(data_AT$group_treatment %in% c("Contestatory", "Collaborative", "Open communication"))
data_AT_treat$group_treatment <- droplevels(data_AT_treat$group_treatment)
tab_xtab(data_AT_treat$final_JR, data_AT_treat$group_treatment, show.col.prc = T, show.summary = F, var.labels = c("Justification rationality", "Treatment"))


# 2) Bivariate relationship between communication formats and Constructive Politics
tab_xtab(data_AT_treat$proposal_final, data_AT_treat$group_treatment, show.col.prc = T, show.summary = F, var.labels = c("Constructive politics", "Treatment"))


# 3) Bivariate relationship between justification rationality and constructive politics
tab_xtab(data_AT$final_JR, data_AT$proposal_final, show.col.prc = T, show.summary = F, var.labels = c("Justification rationality", "Constructive politics"))


# 4) Bivariate relationship between justification rationality and opinion change
tab_xtab(data_AT$opinion_formation, data_AT$final_JR, show.col.prc = T, show.summary = F, var.labels = c("Opinion formation", "Justification rationality"))


# 5) Bivariate relationship between constructive politics and opinion change
tab_xtab(data_AT$opinion_formation, data_AT$proposal_final, show.col.prc = T, show.summary = F, var.labels = c("Opinion formation", "Constructive politics"))








# 6. Regression with weights ----------------------------------------------




# Study 2 Model A ---------------------------------------------------------



# estimate model
at_pol <- multinom(opinion_formation ~ group_treatment + strong_prior + index_corona_knowledge +
                     leftright + trust_gov + trust_exp + index_cognition +
                     index_evaluation + index_accuracy, data = data_AT, weights = data_AT$wght) 

at_pol_N <- nrow(at_pol$residuals) 
at_pol_N # 2110 complete cases

# coef plot study 2 model A
coef_plot_study2_mA <- plot_model(at_pol, type = "est", show.values = T, transform = NULL,
                                  auto.label = F, vline.color = "black", 
                                  value.size = 2.5, dot.size = 1.5, line.size = .5, 
                                  value.offset =.35, order.terms = c(1:24),
                                  axis.labels= c("Accuracy motivation", "Neeed for Evaluation",
                                                 "Need for Cognition", "Trust experts",
                                                 "Trust government", "Leftright",
                                                 "Corona knowledge", "Strong prior",
                                                 "Open Communication", "Collaborative",
                                                 "Contestatory","Information Only"))+
  labs(title = "(A) All Participants", caption = paste0("N = ", at_pol_N))+
  theme_sjplot()+
  scale_y_continuous(limits = c(-3,3))
coef_plot_study2_mA[["layers"]][[1]][["aes_params"]]$size <- 0.3
coef_plot_study2_mA[["layers"]][[1]][["aes_params"]]$alpha <- 0.5

# predicted probabilities
pred_at_pol <- ggemmeans(at_pol, terms = "group_treatment", typical = "mean", condition = c(strong_prior =0, gender_male = 0))
pred_at_pol$x <- factor(pred_at_pol$x, levels = c("Open communication", "Collaborative", "Contestatory", "Information only", "No information (placebo)"))

prob_plot_study2_mA <- pred_at_pol %>%
  ggplot(aes(predicted, x))+
  geom_point()+
  geom_segment(aes(x=conf.low, y =x, xend=conf.high, yend=x))+
  facet_grid(~response.level, scales = "free")+
  xlab("")+
  ylab("")+
  theme_bw()+
  scale_x_continuous(n.breaks = 4, labels = label_percent(accuracy = 1))+
  theme(axis.text = element_text(size =7))


study2_plotA <- coef_plot_study2_mA / prob_plot_study2_mA





# Study 2 Model B ---------------------------------------------------------



# estimate model
at_pol_health <- multinom(opinion_formation ~ group_treatment + strong_prior + index_corona_knowledge +
                            leftright + trust_gov + trust_exp + index_cognition +
                            index_evaluation + index_accuracy, data = data_AT_health, weights = data_AT_health$wght) 

at_pol_health_N <- nrow(at_pol_health$residuals) 
at_pol_health_N # 1319 complete cases

# coef plot study 1 model B
coef_plot_study2_mB <- plot_model(at_pol_health, type = "est", show.values = T, transform = NULL,
                                  auto.label = F, vline.color = "black", 
                                  value.size = 2.5, dot.size = 1.5, line.size = .5, 
                                  value.offset =.35, order.terms = c(1:24),
                                  axis.labels= c("Accuracy motivation", "Neeed for Evaluation",
                                                 "Need for Cognition", "Trust experts",
                                                 "Trust government", "Leftright",
                                                 "Corona knowledge", "Strong prior",
                                                 "Open Communication", "Collaborative",
                                                 "Contestatory","Information Only"))+
  labs(title = "(B) Pro Health", caption = paste0("N = ", at_pol_health_N))+
  theme_sjplot()+
  scale_y_continuous(limits = c(-3,3))+
  theme(axis.text.y = element_blank())
coef_plot_study2_mB[["layers"]][[1]][["aes_params"]]$size <- 0.3
coef_plot_study2_mB[["layers"]][[1]][["aes_params"]]$alpha <- 0.5

# predicted probabilities
pred_at_pol_health <- ggemmeans(at_pol_health, terms = "group_treatment", typical = "mean", condition = c(strong_prior =0, gender_male = 0))
pred_at_pol_health$x <- factor(pred_at_pol_health$x, levels = c("Open communication", "Collaborative", "Contestatory", "Information only", "No information (placebo)"))

prob_plot_study2_mB <- pred_at_pol_health %>%
  ggplot(aes(predicted, x))+
  geom_point()+
  geom_segment(aes(x=conf.low, y =x, xend=conf.high, yend=x))+
  facet_grid(~response.level, scales = "free")+
  xlab("")+
  ylab("")+
  theme_bw()+
  scale_x_continuous(n.breaks = 4, labels = label_percent(accuracy = 1))+
  theme(axis.text.y = element_blank())


study2_plotB <- coef_plot_study2_mB / prob_plot_study2_mB


# Study 2 Model C ---------------------------------------------------------


# estimate model
at_pol_free <- multinom(opinion_formation ~ group_treatment + strong_prior + index_corona_knowledge +
                          leftright + trust_gov + trust_exp + index_cognition +
                          index_evaluation + index_accuracy, data = data_AT_free, weights = data_AT_free$wght) 

at_pol_free_N <- nrow(at_pol_free$residuals) 
at_pol_free_N # 791 complete cases

# coef plot study 1 model C
coef_plot_study2_mC <- plot_model(at_pol_free, type = "est", show.values = T, transform = NULL,
                                  auto.label = F, vline.color = "black", 
                                  value.size = 2.5, dot.size = 1.5, line.size = .5, 
                                  value.offset =.35, order.terms = c(1:24),
                                  axis.labels= c("Accuracy motivation", "Neeed for Evaluation",
                                                 "Need for Cognition", "Trust experts",
                                                 "Trust government", "Leftright",
                                                 "Corona knowledge", "Strong prior",
                                                 "Open Communication", "Collaborative",
                                                 "Contestatory","Information Only"))+
  labs(title = "(C) Pro Civil-Liberties", caption = paste0("N = ", at_pol_free_N))+
  theme_sjplot()+
  scale_y_continuous(limits = c(-3,3))+
  theme(axis.text.y = element_blank())
coef_plot_study2_mC[["layers"]][[1]][["aes_params"]]$size <- 0.3
coef_plot_study2_mC[["layers"]][[1]][["aes_params"]]$alpha <- 0.5

# predicted probabilities
pred_at_pol_free <- ggemmeans(at_pol_free, terms = "group_treatment", typical = "mean", condition = c(strong_prior =0, gender_male = 0))
pred_at_pol_free$x <- factor(pred_at_pol_free$x, levels = c("Open communication", "Collaborative", "Contestatory", "Information only", "No information (placebo)"))

prob_plot_study2_mC <- pred_at_pol_free %>%
  ggplot(aes(predicted, x))+
  geom_point()+
  geom_segment(aes(x=conf.low, y =x, xend=conf.high, yend=x))+
  facet_grid(~response.level, scales = "free")+
  xlab("")+
  ylab("")+
  theme_bw()+
  scale_x_continuous(n.breaks = 4, labels = label_percent(accuracy = 1))+
  theme(axis.text.y = element_blank())


study2_plotC <- coef_plot_study2_mC / prob_plot_study2_mC


# Fig. complete
(coef_plot_study2_mA + coef_plot_study2_mB + coef_plot_study2_mC) /
  (prob_plot_study2_mA + prob_plot_study2_mB + prob_plot_study2_mC)
ggsave("study2_regression_opinion change_logit_wght.png", height = 7, width = 12, dpi = 300)


